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1 . INTRODUCTION 


Under  the  sponsorship  of  the  Office  of  Naval  Research  (ONR) , Nielsen 
Engineering  & Research,  Inc.  (NEAR)  is  conducting  an  ongoing  program  the 
objective  of  which  is  the  testing  of  turbulence  models  using  the  most 
accurate  methods  of  computing  turbulent  flows  now  available.  In  order  to 
limit  the  number  of  possible  sources  of  discrepancy  between  the  predic- 
tions of  a model  and  the  results  of  an  accurate  computation,  the  program 
has  begun  by  looking  at  the  simplest  turbulent  flows  and  the  simplest 
models;  new  features  (geometric  complexity  and/or  physical  phenomena)  will 
be  added  one  at  a time.  This  should  provide  an  expanding  base  of  confi- 
dence to  build  on  and  will,  we  hope,  avoid  some  of  the  difficulties  that 
other  researchers  have  had  in  sorting  out  various  effects  in  turbulent 
flows  and  in  learning  tc  model  them. 

The  basic  approach  to  model  validation  used  is  to  compute  an  accurate 
estimate  of  the  quantity  to  be  modeled  and,  simultaneously,  the  value  that 
the  model  would  predict  for  the  same  quantity.  Comparison,  usually  by 
means  of  a correlation  coefficient,  then  provides  the  information  as  to 
the  validity  of  the  model;  parameters  in  the  model  can  also  be  evaluated. 

This  general  approach  could  be  implemented  in  several  ways . In  the 
first,  exact  Navier-Stokes  simulations  of  turbulent  flows  could  be  used  to 
validate  models  used  in  the  time-average  Navier-Stokes  equations  or  models 
used  in  the  spatially-averaged  equations  solved  in  a large-eddy  simulation. 
In  the  time-average  computation,  the  Reynolds  stresses,  which  are  essenti- 
ally the  long-time  averages  of  products  of  the  fluctuating  components  of 
the  velocity,  need  to  be  modeled.  In  large-eddy  simulation,  on  the  other 
hand,  averages  of  the  products  of  small-scale  components  occur  and  require 
modeling.  By  analogy,  these  are  called  the  subgrid-scale  Reynolds  stresses 
but  they  represent  different  physical  phenomena  than  the  time-average 
Reynolds  stresses.  However,  it  is  believed  that  both  sets  of  quantities 
can  be  modeled  in  similar  ways.  There  is  a difficulty  in  using  exact 
simulation  to  test  time-average  models,  namely  that  only  a few  turbulent 
flows  (all  at  low  Reynolds  number)  are  accessible  to  exact  simulation  and 
this  severely  restricts  what  can  be  done  in  terms  of  time-average  model 
testing. 


In  another  possible  implementation  of  the  general  approach  for  model 
testing,  large-eddy  simulations  could  be  used  to  test  time-average  models; 
a possible  problem  here  is  that  this  testing  is  less  rigorous  because  the 


effects  of  the  subgrid-scale  turbulence  (which  has  to  be  modeled  in  the 
large-eddy  simulation)  are  present  in  the  "baseline"  flow. 

Since  it  is  the  approach  that  provides  the  best  accuracy,  we  have 
chosen  to  begin  by  using  exact  Navier-Stokes  simulations  to  test  subgrid- 
scale  models.  We  believe  that  this  is  the  area  in  which  the  most  informa- 
tion can  be  generated  in  the  shortest  time.  One  can  also  argue  that  the 
same  models  ought  to  be  good  for  both  kinds  of  modeling  (subgrid-scale  and 
time-average)  and  we  are  therefore  indirectly  generating  information  about 
time-average  modeling.  In  later  stages  of  the  program,  we  intend  to  eval- 
uate time-average  models  more  directly  by  comparison  with  large-eddy  simula- 
tion. This  will  be  discussed  further  in  a later  section. 

The  bulk  of  the  work  to  date  has  been  in  the  area  of  the  evaluation 
of  the  validity  of  eddy-viscosity  models  for  the  subgrid-scale  field  in 
homogeneous  isotropic  turbulence  at  relatively  low  Reynolds  numbers.  Some 
study  has  also  been  made  of  the  effects  of  mean  strain  on  this  type  of 
modeling.  The  computer  code  necessary  for  this  work  was  developed  in  the 
first  year  of  this  program  (ref.  1) ; most  of  the  analysis  was  completed  in 
the  second  year  for  which  this  is  the  Technical  Report. 

The  next  section  of  this  report  contains  a brief  description  of  the 
methodology  used  to  do  model  testing.  The  third  section  contains  the  pre- 
sentation and  discussion  of  the  results  generated  to  date,  the  fourth  sec- 
tion contains  the  conclusions  from  this  work,  and  the  last  section  describes 
extensions  of  this  work  to  be  undertaken  in  succeeding  contract  years. 

2 . METHODOLOGY 

As  discussed  above,  our  work  to  date  has  involved  using  the  results 
from  a direct  simulation  of  the  Navier-Stokes  equations  to  test  subgrid- 
scale  models  of  the  eddy-viscosity  type.  The  direct  simulation  used  is 
that  developed  by  R.  S.  Rogallo  (ref.  2)  on  the  ILLIAC  IV  at  the  NASA/Ames 
Research  Center.  This  simulation  calculates  the  evolution  in  time  of  in- 
compressible homogeneous  turbulence  in  the  presence  of  a simple  class  of 
spatially-linear  mean  flows.  A 643  grid  with  spacing  A is  used  which 
allows  direct  simulation  of  flows  with  Taylor-microscale  Reynolds  number 
(Rx)  less  than  about  40;  in  this  range,  the  results  from  the  simulation 
can  be  considered  to  be  exact.  We  are  furnished  a magnetic  tape  contain- 
ing the  exact  velocity  field  u^  on  the  64^  grid  at  an  instant  in  time  as 
calculated  by  this  simulation.  The  tape  is  then  read  into  the  NASA/ARC 
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CDC-7600  computer.  Subsequent  calculations  are  done  on  this  machine  and 
follow  the  approach  developed  by  Clark  (ref.  3) . 

Imagine  that  the  same  homogeneous  incompressible  flow  is  to  be  calcu- 

3 

lated  by  means  of  large-eddy  simulation  on  a coarser  grid  (16  , of  spacing 
Ac  = 4A)  overlaid  on  the  64J  grid.  In  large-eddy  simulation  one  attempts 
to  calculate  the  behavior  of  the  large-scale  components  of  turbulence 
while  the  small  scales  must  be  modeled.  The  large-scale  field  is  obtained 
by  filtering  the  velocity  field;  essentially  this  means  local  averaging 
over  a small  spatial  volume.  The  difference  between  the  instantaneous 
velocity  (u^)  and  the  filtered  value  is  called  the  subgrid-scale  com- 
ponent. That  is, 

ui  = + u|  (1) 

where  u^  is  the  large-scale  (resolvable)  component  defined  by 

u.  = | G(x-x')  ui(x')dx'  (2) 

G is  an  appropriate  filter  function,  the  integration  is  over  the  entire 
flow  field,  and  u!^  is  the  subgrid-scale  field.  The  filter  functions  used 
in  this  work  are  presented  in  Appendix  A.  In  the  equations  describing 
the  large-scale  field  which  are  solved  in  a large-eddy  simulation  (see, 
for  example,  ref.  1)  the  nonlinearity  of  the  Navier-Stokes  equations 
results  in  terms  involving  u^u^  which  play  the  role  of  subgrid-scale 
Reynolds  stresses;  these  are  the  terms  for  which  we  are  evaluating  models. 
Because  we  know  the  exact  value  of  u.  at  each  point  on  the  64 3 grid  from 

* It 

Rogallo's  simulation,  we  can  use  equations  (2)  and  (1)  to  calculate  the 

exact  values  of  u.  and  u! ; we  can  then  calculate  exact  values  for  the 
1 1 3 

subgrid-scale  stress  on  the  64  grid  and  extract  the  appropriate  values 
3 3 

on  the  16  grid.  On  this  16  grid,  we  can  also  calculate  what  various 
eddy-viscosity  type  models  would  predict  for  the  subgrid-scale  stress, 
since  this  type  of  model  represents  this  quantity  as  a functional  of  the 
u^  field.  Note,  however,  that  this  "model  prediction"  is  what  a model 
would  produce  if  it  were  applied  to  the  exact  resolvable  field;  in  an 
actual  large-eddy  simulation,  it  is  applied  to  a resolvable  field  computed 
with  the  use  of  the  model,  and  this  is  somewhat  different.  Our  approach 
does,  nevertheless,  produce  a true  test  of  the  model. 


This  convolution  integral  is  actually  evaluated  using  Fourier  transforms. 
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Since  at  this  point  in  our  procedure  we  possess  exact  and  "modeled" 

3 

values  for  the  subgrid-scale  stress  at  each  point  on  the  16  grid,  we  can 
compare  them  and  get  an  assessment  of  the  accuracy  of  the  models;  correla- 
tion coefficients  are  used  to  obtain  a quantitative  assessment  of  this 
accuracy.  We  can  also  adjust  the  values  of  the  parameters  in  the  models 
to  force  agreement  (in  the  r.m.s.  sense)  with  the  exact  values,  and 
thereby  calculate  what  the  value  of  the  model  parameters  should  be. 

In  the  remainder  of  this  siction,  we  define  some  terms  and  present 
the  models  used.  For  a more  detailed  description  of  the  procedure  briefly 
outlined  above  which  results  in  the  calculation  of  correlation  coeffi- 
cients and  model  parameters,  see  reference  1. 

The  subgrid-scale  stress  is  defined  to  be 

1 


Tij  = 


-s-  6 . . U,  U. 

5 k k 


u,'  u; 


u:u: 
i 1 


(3) 


In  eddy-viscosity  models,  the  subgrid-scale  stress  is  represented  as 
proportional  to  the  rate  of  strain  of  the  resolvable-scale  field 


Tij  ■ 2vT  sij 


where 


S.  . = ^ 
13  2 


3u. 

i 

9x. 


3u . 


3x. 
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The  models  for  vT  evaluated  are: 


1. 


2. 


3. 


2 — — 1/2 

The  Smagorinsky  Model,  v„  = (C  A ) [2S . .S . . ) ' 

1 S d 1]  1] 

filter  length  scale  and  C£  is  the  model  parameter 

2 — — 1/2 

The  Vorticity  Model,  vT  = (CyAa/  (u^u^)  , 

parameter  and  u).  = e.  ..  3u./3x.. 

1 1]K  K J 

The  Kinetic  Energy  Model,  v„  . 

A ^ a 


(4) 


(5) 


where  A=  is  the 

a 


where  Cv  is  the  model 


(C„AJ  (f  u-u-  ) 


1/2 


where  Cg  is 


the  model  parameter  and  1/2  u£u£  is  the  exact  subgrid-scale 
kinetic  energy. 

4.  The  constant-eddy-viscosity  model,  vT  = Cc> 

The  assessment  of  the  accuracy  of  a modeled  quantity  (M)  in  repre- 
senting an  exact  quantity  (X)  is  done  by  means  of  a correlation  coeffi- 
cient 

C(M,X)  = r/T^ — tt*  (6) 


<M2> 


<x2> 
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lated)  to  1 (if  the  model  is  exact  to  within  a multiplicative  constant) . 
These  correlation  coefficients  are  independent  of  the  value  used  for  the 
parameter  in  a given  model.  Values  for  the  parameters  can  be  determined, 
however,  by  forcing  agreement  of  the  root  mean*  square  (r.m.s.)  modeled 
and  exact  values.  For  a further  discussion,  see  Appendix  B. 

While  in  the  foregoing,  the  model  testing  was  discussed  only  in  terms 
of  the  subgrid-scale  stress  itself,  correlation  coefficients  and  model 
parameters  can  actually  be  calculated  on  the  three  levels  used  by  Clark 
(ref.  3) : (1)  the  tensor  level,  where  models  (M^j)  are  compared  directly 

to  the  exact  subgrid-scale  stress  (x^j) ; (2)  the  vector  level,  where 
9Mij/3Xj  is  compared  to  9i^j/3x^;  and  (3)  the  scalar  level,  where  vL 

3Mij/3Xj  is  compared  to  u^  3x^j/3Xj.  Notice  that  at  the  tensor  level,  a 

correlation  coefficient  is  calculated  for  each  stress  component  (6  of 
which  are  independent) ; the  arithmetic  average  of  these  coefficients  gives 
a sense  of  the  correlation  for  the  whole  tensor.  Similarly,  at  the  vector 
level,  there  are  three  correlation  coefficients  (and  their  average) , and 
at  the  scalar  level,  one  correlation  coefficient.  With  respect  to  the 

model  parameters,  at  the  tensor  level  a value  for  the  parameter  can  be 

calculated  for  each  of  the  six  independent  stress  components;  a separate 
value  of  the  parameter  can  be  obtained  for  the  diagonal  and  off-diagonal 
components;  and  still  another  value  of  the  tensor-level  parameter  can  be 
derived  by  applying  the  r.m.s.  criterion  to  the  sume  of  squares  of  the 
stress  components  at  each  grid  point.  Similarly,  three  values  for  the 
parameter  (and  an  "overall"  value)  are  available  at  the  vector  level,  and 
one  at  the  scalar  level. 

All  differentiation  on  the  643  grid  (e.g.,  to  calculate  the  exact 
3 t i n 

value  of  gx  ■*)  is  done  using  Fourier  or  pseudo-spectral  differencing.  The 
differentiation  required  in  the  model  calculations  on  the  163  grid  (e.g., 
to  calculate  S. • or  the  divergence  of  the  model  value  of  the  subgrid-scale 
stress  — g--^-)  is  done  using  any  of  three  methods;  second-order  central 
differencing,  fourth-order  central  differencing,  or  pseudo-spectral  dif- 
ferencing. This  point  will  be  discussed  further  in  the  next  section. 


1 3 

On  the  16  mesh. 
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3.  RESULTS  AND  DISCUSSION 


The  work  in  this  program  is  an  extension  of  that  of  Clark  (ref . 3) . 
He  evaluated  the  eddy-viscosity  models  treated  here  using  one  homogeneous 
isotropic  flow  field  and  one  choice  of  filter.  Our  first  task  was  to 
attempt  to  reproduce  the  calculations  done  by  Clark;  this  is  described  in 
the  next  section.  Following  that,  we  present  various  extensions  of 
Clark's  work  that  we  have  undertaken. 


3.1  Homogeneous  Isotropic  Turbulence 

3.1.1  Comparison  with  the  results  of  Clark.-  Our  calculations  are 
totally  independent  of  Clark's  as  he  did  his  own  direct  simulation  and 
had  his  own  data  reduction  routine.  For  this  comparison  we  did  the  cal- 
culations with  the  same  filter  (a  "box"  filter)  and  filter  width  (Aa  = 8A) 
as  Clark  and  applied  them  in  the  same  way.  However,  some  differences 
remain;  Table  I summarizes  the  parameters  and  methods  used  in  both 
investigations.  There  are  differences  in  detail  in  the  Navier-Stokes 
solvers  used  which  can  result  in  small  differences  in  the  "exact"  velo- 
city fields  produced,  although  the  overall  features  of  both  calculations 
are  in  excellent  agreement  with  each  other  and  with  experiment  as  is  seen 
in  figures  1-3.  In  these  figures,  the  time  dependence  of  the  flow's 
energy,  dissipation  rate  and  Taylor-microscale  Reynolds  number  as  calcu- 
lated by  Rogallo  are  compared  to  values  calculated  by  Clark  and  to  the 
experimental  values  of  reference  4.  The  initial  and  final  values  in 
Clark's  computation  are  shown  as  are  the  values  at  time  step  40;  Clark 
chose  the  velocity  field  at  this  time  step  for  detailed  analysis,  so  we 
have  done  the  same.  The  results  described  below  are  for  the  flow  field 
at  this  time  step. 


The  differences  shown  in  Table  I for  the  calculations  involving  the 
filtered  fields  involve  the  choice  of  coarse  grid  size  and  the  differenc- 
ing method  to  be  used  in  the  model  calculations.  We  chose  a 163  coarse 
grid  because  it  would  result  in  improved  statistics  relative  to  those 
from  the  83  grid  used  by  Clark.  This  is  expected  to  be  a small  effect. 
Also,  we  used  three  differencing  methods  on  our  coarse  grid  whereas  Clark 
used  a single  method  on  his  coarse  grid.  The  effect  on  our  results  of 
changing  the  differencing  method  used  is  discussed  in  detail  in  a later 
section.  In  this  section,  we  compare  Clark's  results  to  our  results 
calculated  using  both  second-order  and  fourth-order  central  differencing; 
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TABLE  I.  SUMMARY  OF 
WORK  AND 

PARAMETERS  AND  METHODS  OF 
THOSE  OF  CLARK  (REF.  3) 

PRESENT 

Clark 

Present  Work 

Navier-Stokes  Solution 

(ref.  2) 

Grid  (Spacing  = A) 

64  3 

643 

Space  Differencing 

Fourth-Order 

Finite  Difference 

Spectral 

Time  Differencing 

Third-Order 

Predictor-Corrector 

Fourth-Ordei 

Runge-Kutta 

Initial  Energy  Spectrum 

From  ref.  4 

Same  as  Clark' 

Initial 

38.1 

Same  as  Clark'; 

Filtered  Fields 

Grid 

83 

163 

Filter 

Box 

Box 

Filtering  Length 

Scale  (A  ) 

d 

8A 

8A 

Model  Derivatives 

Fourth-Order 

Finite  Difference 

Variable 
(see  text) 

11 
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TIME  STEP,  N 


Figure  1.  Turbulence  kinetic  energy  per  unit  mass  as  a 
function  of  time  step.  Note  shifted  origin. 


DISSIPATION  RATE 


TIME  STEP,  N 


Figure  2.  Viscous  dissipation  rate  per  unit  mass  as  a 
function  of  time  step.  Note  shifted  origin. 
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the  truncation  error  for  Clark's  differencing  method  on  his  coarse  grid 
is  bounded  by  the  truncation  errors  for  these  methods  on  our  coarse  grid. 

In  Table  II,  the  correlation  coefficients  and  model  parameters  cal- 
culated in  the  present  work  and  by  Clark  for  the  flow  field  at  time  step 
40  in  the  decay  of  the  homogeneous  isotropic  flow  are  shown.  Examination 
of  this  table  shows  that  the  differences  between  the  calculations  are 
slight.  The  greatest  difference  is  the  lack  (in  the  present  work)  of  an 
increase  in  the  correlation  coefficients  in  going  from  the  tensor  to  the 
vector  level.  This  difference  and  the  small  differences  in  the  parameters 
shown  are  felt  to  be  due  to  small  differences  in  the  velocity  fields  and 
the  different  coarse  grids  used.  In  any  case,  the  differences  are  not 
large  enough  to  modify  Clark's  conclusions,  the  major  one  of  which  is 
that  all  the  models  investigated  perform  essentially  equally  well.  That 
is,  there  seems  to  be  a limit  on  the  accuracy  of  models  of  the  eddy- 
viscosity  type.  An  eddy-viscosity  model  constructed  with  knowledge  of 
the  exact  kinetic  energy  in  the  subgrid-scale  field  predicts  the  energy 
transferred  to  that  field  no  better  than  a constant-eddy-viscosity  model. 
Because  of  this  conclusion,  many  of  the  results  presented  in  the  remain- 
der of  this  report  will  be  for  the  Smagor insky  model  only.  This  model  is 
widely  used  in  large-eddy  simulation  and  the  results  shown  for  it  are 
representative  of  those  derived  for  the  other  models. 

The  comparisons  of  Table  II  have  provided  confidence  in  the  methods 
employed  in  this  study.  We  now  move  on  to  extending  Clark's  work  by 
studying  the  effect  on  the  quality  of  the  models  of  changing  various 
parameters  associated  with  the  flow  fields  studied  and  with  the  modeling 
itself. 

3.1.2  Filter  type.-  In  previous  work  on  large-eddy  simulation  (ref. 
5)  the  choice  of  filter  type  was  found  to  have  only  minor  effects  on  the 
results  of  a simulation.  We  have  investigated  the  Gaussian  filter  as 
well  as  the  "box"  filter  using  our  direct  testing  methods.  Our  calcu- 
lated model  correlation  coefficients  and  parameters  are  shown  for  each  of 
these  filter  types  in  Table  III.  This  comparison  was  done  using  the  same 
flow  field  as  before  (time  step  40  in  Rogallo's  simulation,  R^  = 37)  and 
second-order  central  differences  were  used  in  the  model  calculations. 

The  same  filter  width  (A  = 2A_)  was  used  for  both  filters.  The  results 

a C 

in  Table  III  verify  the  previous  conclusion  that  the  choice  of  filter 
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TABLE  II.  COMPARISON  OF  PRESENT  RESULTS  TO  THOSE  OF  CLARK  (REF.  3) 
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TABLE  III.  EFFECTS 

OF  FILTER 

TYPE  ON  CORRELATION 

COEFFICIENTS 

1 

AND  MODEL 

PARAMETERS 

Average  Correlation 

Model 

IT 

Coef ficient 

Parameter 

. i 

Box 

Gaussian 

Box 

Gaussian 

p 

Tensor  Level 

li 

Smagorinsky 

.33 

.30 

.20 

.18 

r 

Vorticity 

.32 

.30 

.22 

.20 

1 1 

Kinetic  Energy 

.36 

.34 

.16 

.15 

r 

Constant  vT 

.35 

.33 

— 

— 

1. 

f- 

Vector  Level 

1. 

Smagorinsky 

.29 

.27 

.23 

.21 

Vorticity 

.30 

.29 

.25 

.23 

l 

Kinetic  Energy 

.32 

.30 

.21 

.20 

j r 

Constant  vT 

.31 

.30 

— 

1 I ^ 

Scalar  Level 

' 

[j 

Smagorinsky 

.54 

.53 

.17 

1 

.16 

Vorticity 

.55 

.55 

.19 

1 

.18 

1 i; 

Kinetic  Energy 

.58 

.57 

.12 

.11 

u 

Constant  vT 

.56 

.55 

— 

— 

I 


[ 


type  is  a matter  of  the  user's  preference.  Accordingly,  the  results  pre- 
sented in  the  remainder  of  this  report  were  obtained  using  a Gaussian 
filter. 

3.1.3  Filter  width.-  A more  significant  parameter  than  the  filter 
type  is  the  filter  width.  The  filter  width  must,  of  course,  be  large 
enough  to  eliminate  any  structures  that  are  smaller  than  the  computational 
grid.  On  the  other  hand,  if  the  width  is  too  large,  most  of  the  inter- 
esting structure  will  be  eliminated  and  a large  part  of  the  computation 
will  be  wasted.  Therefore,  there  must  be  an  optimum  filter  width. 

Using  the  same  flow  field  as  previously  and  again  using  second-order 
central  differences  in  the  model  calculations,  we  have  investigated  a 
Gaussian  filter  of  varying  width.  Figure  4 shows  the  correlation  coef- 
ficients and  model  parameters  obtained  as  a function  of  filter  width  for 
the  Smagorinsky  model.  This  figure  shows  that  the  correlation  increases 
with  increasing  width  of  the  filter.  The  rate  of  increase  is  rather 
large  at  small  values  of  the  width  but  is  less  at  larger  widths.  It 
thus  appears  that  there  is  a minimum  filter  width  that  should  be  used  in 
large-eddy  simulations.  The  model  parameter  (Cg)  tends  to  be  large  at 
small  filter  widths  and  decreases  as  the  width  is  increased.  There  is 
a range  of  widths  over  which  the  parameter  is  nearly  independent  of  the 
width  and  at  higher  values  of  the  width,  the  parameter  begins  to  increase 
again. 

Some  other  information  can  be  brought  to  bear  on  the  question  of 
filter  width.  It  can  be  shown  (ref.  6)  that  in  homogeneous  turbulence 

£ | HI  dv'  = “2v  J r2  E(k)  dk  (7) 

where  i ( ) dV'  denotes  an  average  over  a spatial  volume  V,  k is  the 

wave  number,  E(k)  is  the  three-dimensional  energy  spectrum  of  the  turbu- 
lence, and  the  right  hand  side  of  (7)  is  e,  the  viscous  dissipation  rate 
of  the  turbulence  kinetic  energy  per  unit  mass.  The  portions  of  this 
dissipation  occurring  in  the  resolvable  and  subgrid-scale  fields  can  be 
assessed  by  making  use  of  E (the  three-dimensional  energy  spectrum  of  the 
resolvable-scale  turbulence) , rewriting  the  right  hand  side  of  (7) 
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Figure  4.  Correlation  coefficient  and  model  parameter  as  a function 
of  filter  width,  Smagorinsky  model,  R^  = 37,  Gaussian  filter. 
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-2v  j k2  E dk  = -2v  j k2  E dk  - 2v  j k2(E-E)dk  (8) 

o o o 

and  taking  the  first  term  on  the  right  hand  side  of  equation  (8)  as  the 
viscous  dissipation  of  the  resolvable  field.  In  figure  5 we  have  plotted 
the  scalar-level  correlation  coefficient  obtained  for  various  filter 
widths  against  the  fraction  of  the  viscous  dissipation  occurring  in  the 
subgrid  scales. 
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Results  are  shown  for  the  flow  field  used  previously  as  well  as  fields 
with  three  lower  values  of  Taylor-microscale  Reynolds  number.  These 
lower-R^  fields  will  be  discussed  more  later.  Although  there  is  some 
scatter  in  this  figure,  it  is  clear  that  the  improved  correlation  as 
A / A increases  shown  in  figure  4 is  accompanied  by  and  is  indi- 
cative  of  an  increase  in  the  fraction  of  dissipation  in  the  subgrid 
scales.  This  correspondence  is  shown  explicitly  in  figure  6 where  the 
fraction  of  dissipation  in  the  subgrid  scales  is  plotted  against  filter 
width  for  the  same  flow  fields.  Consideration  of  figures  4 through  6 
leads  one  towards  recommending  large  values  of  A /A„.  There  is,  however, 
a balancing  factor.  Figure  7 shows  the  fraction  of  the  kinetic  energy  of 
the  turbulence  which  resides  in  the  resolvable  scales  as  a function  of 
filter  width  for  the  same  flow  fields.  This  is  calculated  using  the 
relations 


I |v  C Vj/2,  av 
i [ (u  u /2)  dV 


(9) 


which  follow  from  the  time-average  result  of  reference  7 upon  invoking 
the  ergodic  hypothesis  for  homogeneous  turbulence.  The  ordinate  in 
figure  7 is 
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SCALAR-LEVEL  CORRELATION  COEFFICIENT 


FRACTION  OF  DISSIPATION  IN 
SUBGRID  SCALES 

Figure  5.  Scalar-level  correlation  coefficient  as  a function  of 
the  fraction  of  dissipation  in  the  subgrid  scales. 
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FRACTION  OF  ENERGY 
IN  RESOLVABLE  SCALES 


Aq/Aq 


Figure  7.  Fraction  of  the  turbulence  kinetic  energy  in  the 
resolvable  scales  as  a function  of  filter  width. 


From  figure  7 it  is  clear  that  as  the  filter  width  becomes  large,  very 
little  of  the  energy  in  the  actual  flow  is  actually  captured  in  the 
resolvable  scales  — an  obviously  undesirable  situation. 

Taking  all  of  the  above  factors  into  account  (staying  in  the 
"plateau"  region  of  the  model  parameter,  maximizing  the  model  correlation 
coefficient  while  retaining  a reasonable  portion  of  the  energy  in  the 
resolvable  field) , the  best  filter  width  to  use  in  computing  these  flows 
is  between  two  and  four  times  the  width  of  the  mesh  on  which  the  large- 
eddy  simulation  is  to  be  conducted.  This  range  corresponds,  for  the 
flows  considered  here,  to  requiring  that  at  least  80%  of  the  actual  dis- 
sipation occur  in  the  subgrid-scale  field  and  at  least  10%  of  the  total 
energy  remains  in  the  resolved  field.  These  should  be  useful  guidelines 
for  those  doing  large-eddy  simulations,  but  they  should,  of  course,  be 
checked  in  other  flows. 

3.1.4  Finite  difference  method.-  The  choice  of  filter  type  and  width 
must  obviously  be  made  a priori  by  one  doing  a large-eddy  simulation. 
Another  choice  is  the  differencing  method  to  be  used  in  the  model  calcu- 
lations. In  an  examination  of  the  effect  of  this  choice  on  the  results 
of  a large-eddy  simulation  (ref.  8) , it  was  found  that  the  model  parameter 
obtained  by  matching  the  experimental  energy  decay  rate  in  homogeneous 
isotropic  turbulence  is  essentially  independent  of  the  differencing  scheme 
used. 

We  have  brought  our  direct  testing  methods  to  bear  on  the  question 
of  the  choice  of  differencing  method.  We  calculated  the  model  correlation 
coefficients  and  parameters  using  fourth-order  central  differencing  and 
pseudo-spectral  differencing  and  compared  them  to  the  second-order  central 
differencing  results.  The  correlation  coefficients  are  displayed  in 
Table  IV,  the  parameters  in  Table  V.  In  all  of  these  results,  a Gaussian 
filter  of  width  A /A„  = 2 was  used. 

a C 

In  Table  IV,  we  see  that  the  correlation  coefficient  is  nearly 
independent  of  the  differencing  method  used.  Table  V,  on  the  other  hand, 
shows  that  the  model  parameter  depends  strongly  on  the  differencing 


TABLE  IV. 


EFFECT  OF  DIFFERENCING  METHOD  ON  AVERAGE 
CORRELATION  COEFFICIENTS 


Differencing  Method 


Tensor  Level 

Second 

Order 

Fourth 

Order 

Pseudo- 

Spectral 

Smagorinsky 

. 30 

.30 

.29 

Vorticity 

.30 

.30 

.29 

Kinetic  Energy 

.34 

.33 

.32 

Constant  vT 

.33 

.32 

. 31 

Vector  Level 

Smagorinsky 

.27 

.27 

.24 

Vorticity 

.29 

.29 

.26 

Kinetic  Energy 

.30 

. 30 

.28 

Constant  vT 

.30 

.29 

.27 

Scalar  Level 

Smagorinsky 

.53 

.52 

.50 

Vorticity 

.55 

.54 

.52 

Kinetic  Energy 

.57 

.56 

.56 

Constant  v 

.55 

.55 

.54 

man 


TABLE  V. 


EFFECT  OF  DIFFERENCING  METHOD  ON  MODEL  PARAMETERS 


Differencing  Method 


Second 

Fourth 

Pseudo- 

Tensor  Level 

Order 

Order 

Spectral 

Smagorinsky 

.18 

.16 

.14 

Vorticity 

.20 

.17 

.16 

Kinetic  Energy 

.15 

.13 

.12 

Vector  Level 

Smagorinsky 

.21 

.16 

.13 

Vorticity 

.23 

.18 

.14 

Kinetic  Energy 

.20 

.14 

.10 

Scalar  Level 

Smagorinsky 

.16 

.12 

.10 

Vorticity 

.18 

.14 

.11 

Kinetic  Energy 

.11 

.08 

.06 

method.  In  spite  of  the  results  of  reference  8,  this  is  not  entirely 
unexpected  since  the  differencing  methods  differ  mainly  at  the  high  wave 
numbers  that  provide  a large  part  of  the  strain  in  a turbulent  flow.  The 
difference  between  the  conclusion  reached  here  and  in  reference  8 may  be 
due  to  our  use  of  a directly  simulated  flow  field  to  derive  the  value  of 
the  parameter  whereas  in  reference  8 it  was  determined  using  a flow  field 
calculated  using  the  model.  At  any  rate,  the  results  of  Table  V should 
be  of  interest  as  a guide  to  those  doing  large-eddy  simulation.  Also, 
the  dependence  of  the  model  parameters  displayed  there  is  applicable  to 
the  other  cases  shown  in  this  report  even  though,  in  the  interest  of 
brevity,  we  usually  display  model  parameters  calculated  using  second-order 
differencing  only. 

3.1.5  Reynolds  number . - We  expect  that  the  subgrid-scale  turbulence 
(and  the  model  that  is  used  to  approximate  it)  must  depend  on  a Reynolds 
number.  The  appropriate  velocity  scale  should  be  the  r.m.s.  subgrid-scale 
velocity  fluctuation  and  the  length  scale  should  also  be  one  character- 
istic of  the  subgrid  scale.  The  natural  length  scale  is  the  filter  width 
A which  is  also  used  in  the  models.  The  velocity  scale  is  more  difficult 

a 

to  determine  but  we  expect  that  it  ought  to  be  proportional  to  the  r.m.s. 

strain  rate  S of  the  resolvable  scale  and  the  filter  width.  We  have 

therefore  used  SA./v  = R„„_  as  the  representative  Reynolds  number  of  the 
a Sy  s 

subgrid  scale  and  have  investigated  its  influence  on  the  correlation 

coefficients  and  model  parameters.  The  importance  of  investigating  a 

range  of  R (and  in  particular  extending  the  results  to  low  R ) is 
sy  s sgs 

that  in  a wall-bounded  flow,  Rsgs  will  vary  from  high  to  very  low  values 
near  the  wall.  Because  workers  in  large-eddy  simulation  have  had  dif- 
ficulty in  finding  a subgrid-scale  model  to  be  used  near  a wall,  con- 
siderable importance  is  attached  to  this  investigation. 

The  Reynolds  number  can  be  varied  by  changing  any  of  the  quantities 
contained  in  it.  For  a given  realization  of  a flow,  v is  obviously  fixed 
and,  as  the  filter  width  A is  increased,  the  strain  rate  tends  to  go 

Cl 

down.  When  this  is  coupled  to  the  limited  range  of  A&  for  which  the 
models  are  meaningful  (see  above) , we  find  that  only  a limited  range  of 
Reynolds  number  can  be  obtained  from  a single  realization.  For  this 
reason,  we  asked  Dr.  Rogallo  to  furnish  several  flow  fields  at  different 
energy  levels  in  the  decay  of  homogeneous  isotropic  turbulence  at  dif- 
ferent values  of  R^ . These  variations  give  different  values  of  RsgB*  so 
in  this  way  we  were  able  to  cover  a fairly  wide  range  of  Reynolds  number. 
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Figure  8 shows  the  dependence  on  Ksgs  of  the  scalar-level  correlation 
coefficients  determined  for  the  Smagorinsky  model  using  the  second-order 
central  differencing  method.  The  results  shown  in  this  figure  are  typical 
of  those  obtained  using  the  other  differencing  methods  and  other  models. 

As  just  described,  homogeneous  isotropic  flow  fields  at  several  values  of 
energy  level  and  R^  were  analyzed  with  Gaussian  filters  of  varying  widths, 
subject  to  the  constraints  developed  in  Section  3.1.3:  at  least  80%  of 
the  viscous  dissipation  must  be  in  the  subgrid  scales  and  at  least  10%  of 
the  total  energy  must  remain  in  the  resolved  field.  Although  the  results 
in  figure  8 exhibit  some  scatter,  it  does  appear  that  the  correlation 
coefficient  increases  with  increasing  Reynolds  number.  This  is  expected 
because  the  arguments  used  in  constructing  the  models  are  more  likely  to 
be  valid  at  high  Reynolds  numbers. 

The  dependence  on  Rsgs  of  the  scalar- level  model  parameter  for  the 
Smagorinsky  model  using  second-order  central  differencing  is  shown  in 
figure  9,  for  fourth-order  central  differencing  in  figure  10,  and  for 
pseudo-spectral  differencing  in  figure  11.  The  same  flow  fields  were 
used  in  constructing  these  figures  as  were  used  in  figure  8.  In  these 
figures  it  is  seen  that  the  parameter  in  the  model  displays  a definite 
Reynolds  number  dependence.  Fortunately,  for  each  differencing  method 
used,  it  nearly  reaches  its  high-Reynolds-number  asymptote  within  the 
range  covered.  The  values  of  the  parameter  that  we  find  as  the  high  Rsgs 
asymptotes  are  within  the  range  of  values  that  have  been  used  previously 
(0.1  to  0.24,  refs.  5,  9-11)  to  fit  experiments  at  Reynolds  numbers  above 
those  that  we  were  able  to  cover,  although  agreement  with  Lilly's  theo- 
retical high-R^  value  (0.22,  ref.  12)  is  fairly  poor.  This  last  point 
is  perhaps  indicative  of  the  obvious  dangers  inherent  in  extrapolating 
results  obtained  at  low  Reynolds  numbers  to  the  infinite-Reynolds-number 
limit.  The  computed  values  of  figures  9-11  have  been  fit  by  an  expression 
derived  in  Appendix  C. 

3.1.6  Miscellaneous  investigations.-  In  this  section  we  gather  the 
results  from  brief  investigations  into:  (i)  the  effects  of  deriving  the 
model  correlations  and  parameters  from  different  realizations  of  a tur- 
bulent flow;  (ii)  the  effects  on  the  model  correlations  and  parameters  of 
a modified  form  of  the  Smagorinsky  model  suggested  in  reference  13;  (iii) 
the  effects  of  filtering  on  higher-order  quantities;  and  (iv)  the  use  of 
anisotropic  filters. 
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In  all  of  the  work  described  above,  each  set  of  correlation  coeffi- 
cients and  model  parameters  was  derived  from  a single  realization  of  the 
turbulent  flow  involved.  It  is  obviously  of  interest  to  investigate  the 
effects  on  the  model  correlations  and  parameters  of  using  different 
realizations.  In  the  spectral  simulation  of  the  Navier-Stokes  equations 
of  reference  2,  the  initial  Fourier  velocity  amplitudes  are  selected  with 
random  phases  but  with  the  desired  initial  energy  spectrum  and,  of  course, 
they  satisfy  continuity.  In  this  section,  we  present  the  effects  of 
varying  the  initial  phases  of  the  Fourier  velocity  components . Two  addi- 
tional realizations  of  the  homogeneous  isotropic  turbulent  flow  described 
in  figures  1-3  were  calculated  in  this  way,  and  model  correlation  coeffi- 
cients and  parameters  were  calculated  at  time  step  40  in  each  of  these 
realizations  using  a Gaussian  filter  (A  /A„  = 2) . The  resulting  correla- 
tion  coefficients  calculated  using  second-order  central  differencing  in 
the  model  calculations  are  shown  in  Table  VI,  the  model  parameters  in 
Table  VII.  These  tables  show  that  the  effects  of  initial  conditions  on 
the  correlation  coefficients  and  parameters  are  minimal. 

Based  on  some  studies  using  Burgers'  equation  and  the  Direct  Inter- 
action Approximation,  Love  and  Leslie  (ref.  13)  suggested  modification  of 
the  Smagorinsky  model  by  replacing  the  local  value  by  the  square  root  of 
an  average  of  2S^jS^  over  a large  volume  in  the  calculation  of  the  eddy 
viscosity.  In  their  work,  improved  results  were  obtained  using  this 
non-local  version,  and  they  recommended  testing  this  approach  on  "Navier- 
Stokes  turbulence" . Such  a te^t  was  carried  out  for  a case  of  homogeneous 
isotropic  turbulence  with  R^  = 37  using  a Gaussian  filter  of  width 
Aa/Ac  = 2.  The  extent  of  the  volume  over  which  2S^jS^j  was  averaged  was 
varied.  The  results  in  terms  of  the  scalar- level  correlation  coefficients 
and  model  parameters  for  the  three  different  differencing  schemes  are 
shown  in  figure  12.  In  this  figure,  A^  is  the  filter  width  of  the  "box" 
filter  used  to  do  the  volume  averaging.  It  is  seen  that  the  correlation 
is  only  slightly  improved  by  the  volume  averaging  procedure,  and  that  the 
effect  on  Cs  is  small. 

The  effect  of  filtering  on  higher-order  quantities  was  examined 
using  this  same  flow  field  (R^  = 37)  and  a Gaussian  filter  of  varying 
width.  The  higher-order  quantities  investigated  are  the  resolvable-scale 
velocity-gradient  skewness 
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TABLE  VI . EFFECT  OF  INITIAL  CONDITIONS  ON  AVERAGE  MODEL 
CORRELATION  COEFFICIENTS 

Realization 

Tensor  Level  123 


Smagorinsky 

.30 

.31 

.30 

Vorticity 

.30 

.31 

.30 

Kinetic  Energy 

.34 

.35 

.34 

Constant  vT 

.33 

.33 

.32 

Vector  Level 

Smagorinsky 

.27 

.29 

.26 

Vorticity 

.29 

.30 

.28 

Kinetic  Energy 

.30 

.32 

.29 

Constant  v 

.30 

.31 

.28 

Scalar  Level 

Smagorinsky 

.53 

.56 

.53 

Vorticity 

.55 

.56 

.55 

Kinetic  Energy 

.57 

.59 

.57 

Constant  vT 

.55 

.58 

.55 
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12.  Effect  on  scalar-level  correlation  coefficient 
C of  the  modified  Smagorinsky  model  of  ref.  13. 
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where,  in  this  case,  < > denotes  the  average  on  the  64  J grid  and  the 
subscript  summation  convention  is  suppressed.  In  figures  13  and  14  the 
gradient  skewness  and  flatness,  respectively,  of  the  three  velocity  com- 
ponents are  plotted  as  functions  of  Aa/Ac.  These  results  vary  consid- 
erably from  realization  to  realization,  but  they  show  that  large-eddy 
simulation  cannot  be  reliably  extrapolated  to  zero  filter  width  to  obtain 
skewness  or  flatness.  This  confirms  a result  found  in  reference  5. 

In  wall  bounded  flows,  it  is  necessary  to  use  highly  distorted 
meshes  near  the  boundaries.  In  these  cases,  anisotropic  filters  are 
desirable  and  the  appropriate  length  scale  in  the  subgrid-scale  model  is 
no  longer  obvious.  We  have  taken  an  initial  look  at  the  effect  of  allow- 
ing the  width  of  the  filter  to  be  different  in  each  of  the  three  direc- 
tions. For  this  investigation,  we  used  the  flow  field  with  = 37  and 
chose  the  filter  widths  in  the  three  directions.  We  then  calculated  the 
usual  exact  quantities  (t^,  ^Tij/3xj'  etc.)  using  this  anisotropic  filter 
and  then  examined  the  tensor-level  model  parameters  associated  with  two 

choices  of  length  scale  in  the  Smagorinsky  model.  The  first  length  scale 
1/2 

is  (A  A ) and  the  tensor-level  parameter  associated  with  this  choice 


is  C's 


; ^he  second  is  (A 

ij 


A A ) 
~2  a3 


1/3 


a^  a0  a 


, with  associated  parameter  c" 

sij 


Values  of  C'  and  C"  calculated  using  second-order  differencing  for 
sij  sij 

three  choices  of  anisotropic  filter  are  shown  in  Table  VIII.  Also  shown 
for  comparison  is  the  tensor-level  parameter  CB . associated  with  an 
isotropic  filter  with  A /A_  * 2.  For  each  choice  of  anisotropic  filter 
in  this  table,  shows  considerably  less  variability  than  C'g  , 


indicating  that  (i^  A A ) 
al  a2  a3 


1/3 


is  a better  choice  for  the  length  scale 
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TABLE  VIII.  TENSOR-LEVEL  PARAMETERS  FOR  THE  SMAGORINSKY  MODEL 
ASSOCIATED  WITH  ANISOTROPIC  FILTERS  AND  DIFFERENT 
CHOICES  OF  LENGTH  SCALE 
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when  an  anisotropic  filter  is  used.  For  this  flow  field,  however,  the 
use  of  an  anisotropic  filter  reduces  the  accuracy  of  the  model  from  the 
isotropic  filter  case  as  evidenced  by  the  correlation  coefficients  in 
Table  IX. 


3.2  Homogeneous  Turbulence  in  the  Presence  of  Mean  Strain 

Strain  in  the  mean  field  is  known  to  have  an  important  effect  on 
the  structure  of  the  turbulence.  Presumably,  mean  strain  should  also 
affect  the  structure  of  the  subgrid-scale  turbulence  and  should  affect 
the  modeling  in  some  way.  Since  the  mean  strain  occurs  in  the  largest 
scales,  we  would  expect  that  its  effect  would  be  the  largest  on  the  large 
scales.  This  might  mean,  for  example,  that  in  the  model,  the  mean  strain 
and  the  strain  created  by  che  resolvable  field  should  be  treated  dif- 
ferently. 

The  code  of  reference  2 can  compute  the  case  of  homogeneous  turbu- 
lence acted  on  by  irrotational  mean  strain.  In  this  case,  the  velocity 
field  is  decomposed  as 


ui  = Ui  + ui  + ui  (12) 

where  = ax^,  U2  = bx2,  U3  = cx^  is  the  imposed  field,  and  u^  and  u^ 
are  the  resolvable  and  subgrid-scale  velocity  fields,  respectively.  The 
mean  strain  rates  a,  b and  c are  functions  of  time  only  and  a + b + c = 0 
as  required  by  continuity.  In  the  initial  results  presented  in  this 
section,  the  strain-rate  tensor  which  appears  in  the  eddy-viscosity  models 
is  formulated  as 


1 

2 


(13) 


and  we  have  treated  the  case  of  plane  strain  with  constant  strain  rate, 
i . e . , a=-b=T,  c = 0. 


In  this  kind  of  flow,  two  important  time  scales  are  the  dissipation 
time  scale,  which  is  the  ratio  of  the  turbulence  kinetic  energy  per  unit 
mass  and  the  viscous  dissipation  rate  per  unit  mass 
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< (nL  + uj)  (G.  + uj)  > 
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(14) 


and  the  strain  time  scale, 
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In  our  first  study  involving  this  kind  of  flow,  we  have  investigated 
the  effect  of  the  ratio  of  these  time  scales  on  subgrid-scale  modeling. 
Dr.  Rogallo  ran  his  code  for  three  different  strain  rates  ( r ) , starting 
each  case  with  the  initial  energy  spectrum  used  in  the  isotropic  flows 
discussed  earlier.  The  resulting  values  of  ts/td  are  shown  in  figure  15 
as  a function  of  strain  ratio 


(16) 


This  quantity  is  used  as  sudden  distortion  theory  suggests  that  it  is 
the  controlling  parameter.  The  turbulence  kinetic  energy  per  unit  mass 
and  the  viscous  dissipation  rate  per  unit  mass  for  these  flows  are  shown 
in  figures  16  and  17,  respectively.  The  calculations  were  stopped  at 
roughly  the  same  value  of  strain  ratio,  and  model  evaluation  was  conducted 
using  the  flow  fields  at  the  end  of  the  computation. 

Notice  that  the  strain  ratios  achieved  are  in  the  vicinity  of  1.3 
to  1.4.  This  is  a very  modest  amount  of  straining  (in  fact,  we  see  that 
it  is  insufficient  to  produce  an  increase  in  the  turbulence  kinetic 
energy) , but  it  is  nearly  the  highest  amount  that  can  be  accurately 
simulated  using  the  code  of  reference  2 as  it  is  presently  configured. 

The  reason  is  that  this  program  uses  a coordinate  system  which  deforms 
with  the  mean  strain,  so  the  size  and  shape  of  the  computational  region 
change  with  time.  This  choice  of  coordinate  system  is  necessary  com- 
putationally, but  eventually  the  region  becomes  a "pancake"  and  one  of 
its  dimensions  is  too  small  for  the  computation  to  remain  accurate.  As 
a result,  the  total  strain  is  limited  to  fairly  small  values.  We  hope 
to  investigate  higher  strain  ratios  in  future,  but  this  depends  on  some 
modifications  being  made  to  the  code  of  reference  2.  Higher  mean  strain 
rates  will  also  be  studied. 

The  Smagorinsky-model  correlation  coefficients  and  model  parameters 
calculated  for  the  three  terminal  flow  fields  in  figures  15-17  are  shown 
in  figure  18.  In  these  calculations,  second-order  central  differencing, 
an  isotropic  Gaussian  filter  with  Aa/Ac  = 2,  and  the  rate-of-strain  tensor 
defined  in  (13)  were  used.  The  limiting  values  for  isotropic  turbulence 
(tg/td  -*■  o°)  are  shown  for  comparison.  This  figure  shows  that  for  the 

43 


STRAIN  RATIO 

Figure  15.  Ratio  of  time  scales  as  a function  of  strain  ratio 


Figure  17.  Viscous  dissipation  rate  per  unit  mass  as  a 
function  of  strain  ratio. 
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Figure  18.  Correlation  coefficients  and  model  parameters  for 
the  Smagorinsky  model  as  a function  of  time-scale  ratio. 


range  of  strains  and  time-scale  ratios  covered,  there  is  no  significant 
effect  on  the  subgrid-scale  model.  This  point  is  substantiated  in  Tables 
X (correlation  coefficients)  and  XI  (model  parameters)  where  the  results 
for  the  other  models  are  shown  (again  using  second-order  central  differ- 
encing, an  isotropic  Gaussian  filter  with  Aa/Ac  = 2,  and  equation  (13)). 

The  point  was  made  earlier  that  it  might  be  necessary  for  a subgrid- 
scale  model  to  differentiate  between  the  mean  strain  and  the  resolvable- 
scale  strain.  In  equation  (13) , these  phenomena  are  treated  together. 

It  is  of  interest  to  examine  the  effects  on  the  subgrid-scale  models  of 
using  only  the  subgrid-scale  strain  rate  (equation  (5))  where  a strain- 
rate  tensor  is  required  in  the  models.  For  the  highest  mean  strain  rate 
treated  thus  far,  however,  the  r.m.s.  strain  rate  using  equation  (5)  is 
only  12%  less  than  that  calculated  using  equation  (13)  , so  only  small 
effects  should  be  discernible.  In  figure  19,  the  Smagorinsky  scalar-level 
correlation  coefficients  and  model  parameters  calculated  using  equation 
(5)  (solid  symbols)  are  compared  for  the  different  time-scale  ratios  to 
those  from  figure  18  which  were  calculated  using  equation  (13)  (open 
symbols) . As  expected  for  the  range  of  mean  strain  rates  treated,  the 
effect  of  using  equation  (5)  instead  of  equation  (13)  is  negligibly  small. 
However,  the  divergence  between  the  solid  and  open  symbols  increases  as 
tg/t^  decreases,  so  this  question  should  be  re-investigated  when  larger 
mean  strain  rates  are  studied. 

A final  investigation  was  made  using  one  of  these  strained  flow 

1/3 

fields.  In  an  earlier  section,  the  choice  of  (A,  A A_  ) ' was  found  to 

al  a2  a3 

be  the  appropriate  length  scale  for  use  with  an  anisotropic  filter.  This 
conclusion  was  derived  from  isotropic  turbulence,  however,  and  its  use 
was  found  to  result  in  slightly  lower  correlation  coefficients  than  were 
calculated  using  an  isotropic  filter.  This  conclusion  was  reevaluated 
in  the  flow  field  with  tg/td  = 1.18  by  comparing  the  results  from  figure 
18  (where  an  isotropic  Gaussian  filter  with  A /A  = 2 was  used)  with 
those  calculated  using  the  anisotropic  filter  described  below  (equation 
(13)  was  used  in  both  cases) . Because  the  mean  strain  ratio  in  each 
coordinate  direction  is  known  for  this  flow  field,  we  can  choose  the 
filtering  length  scales  such  that  they  remain  a constant  multiple  of  the 
length  of  the  sides  of  a strained  fluid  element.  That  is. 
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ts/td 

1.18 

2.10 

3.34 

00 

Strain  Ratio 

1.32 

1.33 

1.42 

1.00 

1 

Tensor  Level 

1 

Smagorinsky 

.31 

.29 

.29 

.30 

1 

1 

Vorticity 

.30 

.29 

.27 

.30 

; 

Kinetic  Energy 

.34 

.33 

.32 

.34 

[1 

Constant  vT 

.32 

.31 

.30 

.33 

1 

> 

Vector  Level 

I 

1 

1 

Smagorinsky 

.27 

.24 

.23 

.27 

1 

Vorticity 

.29 

.27 

.24 

.29 

• 

1 

Kinetic  Energy 

.30 

.28 

.27 

.30 

i 

« 

- 

Constant  vT 

Scalar  Level 

.29 

.26 

.25 

.30 

j ! 
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Smagorinsky 

.49 

.47 

.45 

.53 
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Vorticity 

.52 

.52 

.48 

.55 

Kinetic  Energy 

.54 

.53 

.51 

.57 
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Constant  vT 
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.51 

.48 

.55 
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TABLE  XI.  EFFECT  OF  IRROTATIONAL  PLANE  STRAIN  ON  MODEL  PARAMETERS 


ts/td 

1.18 

2.10 

3.34 

OO 

Strain  Ratio 

1.32 

1.33 

1.42 

1.00 

Tensor  Level 

Smagorinsky 

.15 

.17 

.17 

.18 

Vorticity 

.18 

.20 

.20 

.20 

Kinetic  Energy 

.13 

.15 

.16 

.15 

Vector  Level 

Smagorinsky 

.18 

.20 

.19 

.21 

Vorticity 

.22 

.22 

.22 

.23 

Kinetic  Energy 

.18 

.19 

.19 

.20 

Scalar  Level 

Smagorinsky 

.14 

.15 

.15 

.16 

Vorticity 

.16 

.17 

.17 

.18 

Kinetic  Energy 

.10 

.11 

..  -11 

.11 

□ MEAN  STRAIN  RATE  INCLUDED, 
EQUATION  (13) 

■ RESOLVABLE- SCALE  STRAIN 
RATE  ONLY,  EQUATION  (5) 
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Figure  19.  The  effect  of  including  the  mean  strain  in  the 
strain- rate  tensor  used  in  the  Smagorinsky  model. 
Scalar-level  results  are  shown. 
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The  correlation 


Under  these  conditions,  notice  that  (A  = A A 

1 a2  “3 

coefficients  and  tensor- level  model  parameters  for  the  Smagorinsky  model 
calculated  using  second-order  differencing  for  this  investigation  are 

associated 


ij 


shown  in  Table  XII.  The  tensor-level  model  parameter  C" 

1/3 

with  (A  A A ) and  the  anisotropic  filter  defined  by  equation  (17) 
al  a2  a3 

is  seen  to  have  less  variability  than  c'  , again  demonstrating  that 
1/3  *3 

(A  A A ) is  the  appropriate  length  scale  for  use  with  an  anisotropic 
al  a2  a3 

filter.  In  contrast  to  the  results  derived  using  isotropic  turbulence 

however,  C"  also  has  less  variability  than  C_  (associated  with  the 
s . . s.. 

IJ  13 

isotropic  filter) , and  the  anisotropic-filter  correlation  coefficients 

are  slightly  improved  over  those  for  the  isotropic  filter.  While  these 

1/3 

results  demonstrate  that  (A,  A A ) is  the  length  scale  of  choice  in 

1 a2  a3 

flow  fields  with  small  anisotropy,  we  must,  of  course,  reserve  final 
judgment  on  this  point  until  we  have  studied  fields  with  higher  strain. 

As  previously  mentioned,  this  is  planned  pending  modifications  to  the 
code  of  reference  2. 


4.  CONCLUSIONS 


The  studies  we  have  made  of  eddy-viscosity  type  subgrid-scale  models 
as  applied  to  homogeneous  isotropic  turbulence  and  to  homogeneous  turbu- 
lence in  the  presence  of  irrotational  plane  mean  strain  have  led  to  the 
following  conclusions: 

1.  The  earlier  analysis  of  reference  3 of  eddy-viscosity  models 

using  a box  filter  and  a single  homogeneous  isotropic  flow  field 
has  been  repeated  with  only  small  differences  in  the  results. 
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TABLE  XII.  CORRELATION  COEFFICIENTS  AND  PARAMETERS  FOR  THE 
SMAGORINSKY  MODEL  FOR  A FLOW  FIELD  WITH  I RROTATIONAL 
PLANE  STRAIN.  EFFECTS  OF  AN  ANISOTROPIC  FILTER 
AND  DIFFERENT  LENGTH  SCALES 


Correlation  Coefficients 


Tensor  Level 

i*  j 

Isotropic 

Filter 

. 

Anisotropic 

Filter 

1.1 

.47 

.54 

2,2 

.32 

.46 

3,3 

.27 

.33 

1,2 

.31 

.32 

1,3 

.32 

.28 

2,3 

.23 

.33 

Average 

.31 

.35 

Vector  Level 

i 

1 

.40 

.39 

2 

.10 

.13 

3 

.30 

.35 

Average 

.27 

.29 

Scalar  Level 

Model  Parameters 

Isotropic  Anisotropic 
Filter  Filter 


8ij 

Cs.. 

13 

c;.. 

13 

.15 

.11 

.15 

.13 

.18 

.14 

.16 

.15 

.15 

.17 

.17 

.17 

.19 

.16 

.18 

.16 

.19 

.16 
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The  major  conclusion  of  reference  3 is  verified:  all  of  the 
eddy  viscosity  models  tested  perform  about  the  same  and  demon- 
strate a modest  level  of  correlation  with  the  exact  results. 


' 

I 

| 

I 


. 


2.  The  choice  of  filter  type  has  only  a minor  influence  on  the 
accuracy  of  an  eddy-viscosity  subgrid-scale  model.  The  choice 

of  filter  type  to  be  used  in  a large-eddy  simulation  is  therefore 
a matter  of  the  user's  preference. 

3.  There  appears  to  be  an  optimum  value  for  filter  width.  For  the 
flows  studied,  this  width  is  from  two  to  four  times  the  spacing 
of  the  grid  on  which  the  large-eddy  simulation  is  to  be  con- 
ducted. Outside  this  range,  the  accuracy  of  the  model  deteri- 
orates badly  or  the  simulation  contains  too  small  a proportion 
of  the  flow's  energy  to  be  meaningful.  This  also  verifies  con- 
clusions reached  earlier  by  others. 

4.  The  choice  of  differencing  method  used  in  the  model  calculations 
has  essentially  no  effect  on  the  accuracy  of  the  models  as 
assessed  by  the  correlation  coefficients  we  calculated.  The 
model  parameters,  on  the  other  hand,  are  influenced  in  an 
important  way  by  the  differencing  method  used.  This  conclusion 
is  at  variance  with  what  others  have  found. 

_ 2 

5.  The  subgrid-scale  Reynolds  number  R = SA  /v  seems  to  charac- 

Sg  S cl 

terize  the  subgrid  scale  in  the  flows  considered.  Also,  the 
subgrid-scale  models  investigated  are  more  accurate  at  high 
values  of  subgrid-scale  Reynolds  number  than  at  low  values. 
Because  of  the  lack  of  a more  accurate  alternative,  however, 
their  use  in  the  low-Reynolds-number  range  is  required  near 
the  wall  in  wall-bounded  flows.  Towards  this  end,  a relation 
has  been  found  that  represents  the  variation  with  subgrid-scale 
Reynolds  number  of  the  scalar-level  parameter  of  the  Smagorinsky 
model . 

6.  The  effects  of  deriving  model  correlation  coefficients  and 
parameters  from  different  realizations  of  a turbulent  flow  are 
negligible. 

7.  The  modified  form  of  the  Smagorinsky  model  suggested  in  refer- 
ence 13  involving  the  use  of  an  eddy  viscosity  averaged  over  a 


A.  » 
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large  volume  is  not  found  to  result  in  significantly  improved 
accuracy  and  the  effect  on  the  model  parameter  is  small. 


8.  Higher-order  statistical  quantities  (such  as  velocity-derivative 
skewness  or  flatness)  for  the  resolvable-scale  field  in  a 
large-eddy  simulation  cannot  be  reliably  extrapolated  to  zero 
filter  width  to  obtain  the  unfiltered  values. 


9.  For  the  modest  strain  rates  and  strain  ratios  investigated  in 

this  study,  no  appreciable  effects  of  mean  strain  on  the  subgrid- 
scale  models  are  apparent. 


10.  In  cases  where  it  is  desirable  to  use  an  anisotropic  filter,  the 

1/3 

appropriate  length  scale  appears  to  be  (Aa  Aa^Aa^)  . This 
conclusion  is  based  on  flows  with  small  anisotropy  and  must  be 
reexamined  in  flows  with  large  strain. 


5.  FUTURE  DIRECTIONS 

The  use  of  eddy-viscosity  models  in  strained  homogeneous  turbulence 
will  continue  to  be  investigated.  It  is  anticipated  that  modifications 
to  the  code  of  reference  2 will  be  made  that  will  allow  calculation  of 
flows  with  larger  strain;  when  this  is  done,  effects  of  smaller  tg/td 
and  larger  strain  ratios  will  be  investigated  and  the  use  of  anisotropic 
filtering  explored  in  more  depth. 

Another  important  effort  to  be  pursued  in  the  future  is  the  investi- 
gation of  the  validity  of  models  which  use  differential  equations  for 
the  subgrid-scale  stress,  i.e.,  full  Reynolds  stress  models.  Most  of  the 
terms  involved  cannot  be  measured  experimentally,  so  the  model  testing 
we  will  do  should  provide  information  of  considerable  interest  to  those 
doing  large-eddy  simulation.  The  model  testing  will  proceed  along  the 
same  lines  as  is  described  in  this  report,  but  the  calculation  of  many 
new  terms  is  involved.  The  computer  program  necessary  to  do  these 
calculations  is  being  written  at  the  present  time. 

On  a longer  time  scale,  we  intend  to  begin  a new  direction  for  this 
program:  the  testing  of  time-average  models  using  large-eddy  simulation 

as  the  baseline.  The  underlying  premise  of  this  approach  is  that  large- 
eddy  simulation  is  a sufficiently  accurate  representation  of  the  actual 
flow  in  question  to  allow  its  use  as  the  standard  against  which  the 
(time-average)  models  are  to  be  compared.  This  premise  is,  of  course. 


not  known  to  be  true  in  general  at  this  point.  However,  under  certain 
circumstances,  the  premise  is  reasonable.  For  instance,  if  the  subgrid- 
scale  stresses  are  small  compared  to  the  resolvable-scale  stresses  in  a 
particular  simulation,  the  accuracy  of  the  subgrid-scale  stress  model  used 
is  not  critical  to  the  accuracy  of  the  large-eddy  simulation,  and  the 
resolvable-scale  field  computed  in  this  circumstance  can  be  assumed  to 
be  an  accurate  representation  of  an  actual  flow  field.  This  situation 
could  be  abetted  by  lowering  the  Reynolds  number  of  the  simulated  flow, 
for,  as  the  Reynolds  number  decreases,  the  range  of  scales  in  the  flow 
decreases  and  the  ability  to  "capture"  them  with  a given  grid  is  in- 
creased. Put  in  another  way,  as  the  Reynolds  number  decreases  for  a given 
computational  grid,  the  demands  on  the  subgrid-scale  model  are  decreased 
and  the  large-eddy  simulation  approaches  the  accuracy  of  a direct  simula- 
tion. Therefore,  for  some  flows  in  some  Reynolds  number  range,  large- 
eddy  simulation  will  be  an  acceptable  standard  for  model  testing. 
Furthermore,  this  Reynolds  number  range  will  be  above  that  for  which  a 
direct  simulation  is  possible.  The  precise  directions  of  this  work  and 
the  flows  for  which  model  evaluation  will  be  conducted  depend  on  future 
developments,  so  no  further  details  will  be  presented  here. 


APPENDIX  A 


DEFINITION  OF  THE  FILTER  FUNCTIONS  USED  IN  THIS  WORK 


Two  filter  functions  are  used  in  this  work.  The  first  is  the  "box" 
or  "top  hat"  filter 

3 


|x.  - x'|  < Aa/2,  i = 1,2,: 


G (x  - x')  = 


(A.  1) 


0 otherwise 

The  other  is  the  Gaussian  filter 


<*i  - *i>2 


g (x  - x')  = n 
i=l 


/y/iT 


a . 

l 


(A. 2) 


where  Y = 6 and  the  generalized  (anisotropic)  form  is  shown. 


Fourier  transform  methods  were  used  to  evaluate  the  convolution 
integrals  necessary  to  calculate  filtered  quantities.  In  these  pro- 
cedures, the  discrete  transform  of  the  box  filter  as  calculated  by  a 
system-provided  FFT  was  used,  whereas  the  continuous  transform  of  the 
Gaussian  filter  was  used.  These  choices  were  made  for  convenience  and 
are  not  felt  to  exert  any  influence  on  the  results  obtained. 
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APPENDIX  B 


STATISTICAL  MODEL  TESTING 


The  kind  of  model  testing  we  are  doing  in  this  report  is  one  of  the 
simplest  types;  we  are  checking  the  proportionality  of  two  quantities. 

We  denote  the  exact  stress  by  t and  the  model  as  m = ca  where  c repre- 
sents the  free  parameter  and  o the  function  to  be  tested.  (The  component 
indices  play  no  important  role  here  and  will  be  dropped.)  The  question 
is  how  accurately  the  model  represents  actuality.  The  usual  measure  of 
difference  between  the  model  and  the  exact  result  is  the  mean  square 
deviation 


1 N 

- y 

M L 

N i=l 


cck) 


(B.l) 


where  the  index  i represents  a given  trial  and  N is  the  sample  size. 
In  our  c?*se  N ■ 262,144  but  it  is  difficult  to  argue  that  these  are 
independent  samples  because  there  is  obviously  spatial  correlation  in 
the  data. 

.2 


Minimization  of  the  residual  6' 

T t .a . 
c = i 

1 l a? 


gives 


(B.  2) 


The  reason  for  the  subscript  on  c will  be  obvious  below, 
coefficient  used  as  a measure  of  accuracy  is 

ZTi°i 


P = 


In  terms  of  these,  the  residual 

,2 


1 K) 

<52  is 


172 


The  correlation 


(B.  3) 


6 


(1  - 


P2> 


v 2 
N 


(B.  4) 


which  shows  that  the  root-mean- square  deviation  is  large  even  for  the 
largest  correlations  found  in  this  report. 

Note  that  the  parameter  defined  by  Eq.  (B.2)  differs  from  the  one 
we  have  used  in  the  text.  Equation  (B.2)  gives  the  parameter  which 
produces  the  minimum  deviation  6,  i.e.,  the  best  fit  of  the  model  to 
the  data  in  the  mean-square  sense.  However,  this  parameter  varies  more 
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n 


WdfflSMiBfittE 


with  changing  filter  width  (or  subgrid-scale  Reynolds  number)  than  does 
the  paramet  obtained  by  matching  the  mean-square  values  of  the  model 
and  actual  values: 


c2  = 


v 2 


L i 


1/2 


Cj/P 


(B.5) 


This  gives  a mean  square  deviation 


2(1  - p) 


r 2 

K 


N 


2 

L+p 


(B.  6) 


This  method  of  evaluating  the  parameter 


which  is  always  greater  than  6 
places  heavier  weight  on  the  large  values  of  t and  a than  does  the  first 
method . 


It  appears  that  the  critical  issue  in  constructing  a model  is  that 
the  dissipation  of  kinetic  energy  be  large  in  the  model  where  the  actual 
dissipation  is  large  and  vice  versa.  The  choice  that  leads  to  c2  seems 
to  do  this  better  than  the  first  choice,  gives  a value  in  better  agree- 
ment with  the  one  determined  semi-empirically , and  has  less  variation 
with  RggS-  It  has  therefore  been  adopted  here.  We  will  look  into 
statistical  methods  of  model  testing  more  deeply  in  the  future. 


APPENDIX  C 


REYNOLDS  NUMBER  EFFECTS  ON  MODELS  FOR 
SUBGRID-SCALE  (SGS)  TURBULENCE 

An  important  parameter  in  any  flow,  the  Reynolds  number  can  be 

expected  to  have  an  important  effect  on  the  nature  of  the  subgrid-scale 

turbulence.  It  is  important  to  note  that  the  most  important  SGS  eddies 

are  the  largest  ones,  i.e.,  the  ones  that  are  just  smaller  than  the 

resolvable  length  scale.  Thus,  the  important  length  scale  is  likely  to 

be  the  filter  size  A&.  The  velocity  scale  q for  the  SGS  can  be  expected 

to  be  proportional  to  the  strain  S in  the  large-eddy  field  since  SGS 

production  is  proportional  to  it.  Thus  we  expect  that  q ~ SA  and  the 

— 2 a 

proper  Reynolds  number  for  the  SGS  is  R = SA  /v. 

sgs  d 

At  high  Reynolds  numbers,  the  largest  SGS  eddies  receive  energy 
from  the  resolvable  eddies  and  lose  it  to  still  smaller  ones  which  are 
viscously  dissipated.  From  the  point  of  view  of  the  largest  SGS  eddies 
the  transfer  to  the  smaller  scales  is  a dissipation  or  energy  drain.  If 
we  assume  that  the  production  and  dissipation  rates  of  these  eddies  are 
in  balance  at  every  point  in  space  and  time  (a  reasonably  good  approxima- 
tion in  the  homogeneous  flows  we  consider  in  this  report)  then  it  is 
possible  to  show  that  the  Smagorinsky  model  is  a good  approximation  at 
high  SGS  Reynolds  numbers  and  the  parameter  C should  be  independent  of 

■W 

At  low  Reynolds  numbers  viscosity  becomes  important  for  the  largest 
SGS  eddies  and,  perhaps,  even  for  the  resolvable  ones.  Under  these 
circumstances,  it  is  more  difficult  to  find  a good  model  for  the  SGS 
but  we  can  adopt  the  Smagorinsky  model  as  one  which  has  many  of  the 
properties  we  seek  including  symmetry  and  dissipative  behavior.  However, 
we  expect  that  the  parameter  in  the  model  will  become  a function  of 
Reynolds  number.  The  effect  is  similar  to  the  viscous  damping  that 
occurs  near  a solid  surface.  The  effect  of  viscosity  is  to  reduce  the 
strength  of  the  largest  SGS  eddies  (the  transfer  mechanism  that  dominates 
at  high  RBg8  still  operates  as  well) . In  fact,  Rsgs  + 0 is  equivalent  to 
q + 0 in  which-  case  the  SGS  should  have  no  effect  on  the  resolvable  field 
Thus  Cfl  -*■  0 as  Rsg8  "*■  0 and  linear  proportionality  is  a reasonable 
assumption. 


I 

I 

I 


A functional  form  which  has  the  proper  behavior  at  both  high  and 


and  we  have  used  this  to  fit  our  results.  A "derivation"  of  this  rela- 
tionship based  on  the  equations  is  possible  but  no  more  enlightening 
than  the  heuristic  arguments  we  have  given.  The  parameter  a is  essenti- 
ally a Reynolds  number  at  which  viscous  effects  become  important  (similar 
to  the  viscous  damping  parameter  A+  of  boundary  layer  modeling)  and  the 
value  of  ~20  found  by  curve  fitting  seems  reasonable. 
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